Enhanced Diffusion of a Needle in a Planar Course of Point Obstacles 
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The transport of an infinitely thin, hard rod in a random, dense array of point obstacles is investigated by 
molecular dynamics simulations. Our model mimics the sterically hindered dynamics in dense needle liquids. 
Transport becomes increasingly fast at higher densities, and we observe a power-law divergence of the diffusion 
coefficient with exponent 0.8. This phenomenon is connected with a new divergent time scale, reflected in a 
zigzag motion of the needle, a two-step decay of the velocity-autocorrelation function, and a negative plateau in 
the non-Gaussian parameter. Finally, we provide a heuristic scaling argument for the new exponent. 

PACS numbers: 05.20.Jj, 51.10.+y, 45.05.-x 



Model liquids of infinitely thin, hard needles have triv- 
ial equilibrium properties, but exhibit complex dynam- 
ics. Their equilibrium statistical mechanics is trivial, 
since the excluded volume is zero and all configurations 
are permitted and equally likely; in particular, the equa- 
tion of state and the pair distribution function coincide 
with those of an ideal gas. But the mutual steric hin- 
drance of moving needles gives rise to strong dynamic 
correlations for dense systems. The dynamically in- 
duced steric interaction can even generate dynamic ar- 
rest: a glass transition of needles has been observed in 
simulations |[l|j_0] and is rationalized within a micro- 
scopic theory pi]. Due to the absence of static correla- 
tions, this glass transition is of very different origin than 
in structural glass formers. 

The close relationship between liquids and granular 
gases 1 4] suggests similar phenomena also for granular 
needles and fibers, which constitute highly anisotropic 
particles with convex shape. There, interesting slow 
relaxation phenomena with algebraic tails were pre- 
dicted l^, and are likely to be verified soon as experi- 
mental progress continues |7, 8]. 

In early molecular dynamics simulations of a needle 
liquid, a surprising behavior of the center-of-mass (CM) 
diffusion was found by Frenkel and Maguire |0], who 
conjectured a divergence of the self-diffusion coefficient 
Dcm with increasing density n. In the dilute regime, 
the diffusion coefficient decreases as 1 /n in accord with 
Enskog theory, but it starts growing beyond a crossover 
density. In a theoretical attempt to explain this observa- 
tion, they predicted D^m ~ An increase of D^m by 
a factor 2 above its minimum value was confirmed a few 
years later [10] and also in a more recent simulation flU]. 
Such a behavior contradicts the experience that transport 
becomes slow in dense, complex liquids. 

In this Letter, we present extensive molecular dynam- 
ics simulations for a single needle exploring a planar 
course of impenetrable, frozen point obstacles. This 
model again possesses trivial equilibrium properties, 
but exhibits intriguingly complex dynamics as we shall 
demonstrate below by a surprising enhancement of CM 
diffusion. Since the needle is confined by the obsta- 



cles to a narrow tube, its motion is essentially one- 
dimensional leading to a strong suppression of rotational 
diffusion [12]. The same underlying physics occurs also 
in dense needle liquids, where the tube is formed dynam- 
ically by the surrounding needles. 

The needle is characterized by its length L, mass m, 
and moment of inertia /; we use a homogeneous mass 
distribution, / = mL^/12. The degrees of freedom 
(DOFs) of the needle encompass the CM position R and 
the unit vector of orientation u, parametrized by a sin- 
gle angle in our planar model; the conjugate momenta 
are mv and luj, see Fig. [T] We choose the collisions to 
be elastic and frictionless; then the total kinetic energy is 
conserved and its value merely sets the overall time scale 
of the problem. We distribute the obstacles randomly 
and independently with average number density n; then 
the only dimensionless control parameter of the model is 
the reduced density, n* := nL?. A second length scale 
^ := n^^/^ is given by the typical distance between ob- 
stacles; the crossover from the dilute to the constrained 
regime is expected when ^ competes with L. 

We attach an instantaneous, body-fixed, right-handed, 
ortho-normalized reference frame (u, u^) to the nee- 
dle; it is convenient to define the two-dimensional vec- 
tor product a X b := det(a, b), hence u x = 1. 
The laws of collision for a needle hitting a point obstacle 
follows from two assumptions: (i) the collision is elas- 
tic and (ii) the momentum transfer occurs perpendicu- 
larly to the needle, Ap = Apu^. Simultaneously, a 
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FIG. 1 : Geometry and notation for the needle hitting an obstacle. The 
momentum transfer Ap || ux is uniquely determined for a friction- 
less, elastic collision; it induces also a torque changing the angular 
velocity of the needle. 
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FIG. 2: Simulation results for the diffusion coefficients of the needle. 
Solid lines show the low-density behavior, D ~ 1 /n, dashed lines 
are fits to the asymptotic high-density behavior, and the dotted line 
indicates slope 1/2 for comparison with a scaling prediction. 



collision generates a torque on the needle proportional 
to the lever arm a, given by the distance between the 
collision point and the axis of rotation, which yields a 
total change of angular momentum au x Ap = aAp. 
The post-coUisional velocities Vf = Vi + Ap/m and 
u}f = cji + aAp// follow from the pre-collisional ones, 
Vi and cjj, imposing conservation of kinetic energy. 



Ap = 2m/ 



Vj X u — LOia 
I + ma^ 



(1) 



Our simulations of the molecular dynamics are based 
on an event-driven algorithm |[l3, [3]; the challeng- 
ing part is collision detection for objects of zero width. 
This is solved in an efficient and robust way by apply- 
ing the interval Newton method for root finding [14]. 
It determines the zeros of the time-dependent distance, 
d{t) = [R(t) — rj] X u(t), between the needle and a par- 
ticular obstacle at Vj. For each density, we have run at 
least 300 trajectories, each covering up to 10^ collisions 
in the dense regime. 

The collisions conserve the kinetic energy of the nee- 
dle and redistribute it continually and equally between 
the two translational and one rotational DOF. The phase 
space trajectories are chaotic as in the case of the Lorentz 
model [15, 16], and we have numeric evidence that the 
system is ergodic at all investigated densities; in particu- 
lar, equipartition holds, m(w'^^ /2 = I (uj'^). Further- 
more, the mean collision rate follows the low-density 
result, T~^^ oc n*, also at high density; we have mea- 
sured the prefactor to t^^^ = 0.845 n*w/L, where v de- 
notes the root mean-square velocity. In the low-density 
regime, the rotational and translational diffusion coef- 
ficients obey Dmi^Dcm ~ "^0011 as expected from the 
molecular chaos hypothesis. 



The dense array of obstacles leads to the formation 
of an effective tube caging the needle. Following the 
conventional paradigm, transport should be drastically 
suppressed. Indeed, the rotational diffusion coefficient 
obeys the predicted Doi-Edwards scaling [17] at high 
densities, D^^x ~ (^*)~^ [12]. While the rotational 
diffusion slows down with increasing obstacle density, 
the translational diffusion D^m exhibits a minimum at 
n* f« 5 followed by a rapid increase, see Fig. [2l Within 
the density range covered by our simulations, Dcm is en- 
hanced by more than two decades. Our data suggest a 
power-law divergence. 



in 



for n* 



00, 



(2) 



with an exponent of ^ = 0.80 ± 0.05. In a recent 
simulation for a planar "shish kebab" model with soft 
spheres [18], also a slight increase of Dcm has been ob- 
served. The enhancement of diffusion, however, is cut 
off at moderate densities by the competition with the 
slow dynamics induced by the excluded volume [19]. 

A typical trajectory exhibits a zigzag pattern and is 
mainly composed of long, slightly bent segments, see 
Fig. [3] and the Supplementary Movie [20]; the motion of 
the needle carves out narrow channels. This straight mo- 
tion is occasionally interrupted and the needle rests for a 
short while, until it either continues its way or reverses 
its direction. The latter is reflected in a change of sign 
of the longitudinal velocity v^\, where we split the CM 
velocity into longitudinal and transverse components in 
the co-moving reference frame, v = W||U + v^u^. Such 
a velocity reversal is accompanied by sharp cusps in the 
trajectory, and the general form of the trajectories is rem- 
iniscent of the non-holonomic dynamics of a skate on 
ice jHH] with its local constraints due to friction. 

The time series of the kinetic DOFs in Fig. |4] exempli- 
fies that the transversal and rotational momenta fluctuate 
on the scale of individual collisions, while the longitudi- 
nal velocity fluctuates on a much slower time scale Tz^. 
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FIG. 3: Snapshots of the needle at n* — 120. The time interval 
between subsequent exposures. At = 1.8 L/v, covers about 180 
collisions; colors encode the longitudinal velocity, v«/vmiix- 
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FIG. 4: Time series of the kinetic degrees of freedom for n* = 120. 
Sections of a trajectory may be divided into different states: "run- 
ning" (A), "rattling" (B) and "running backwards" (A), see lower 
panels. The right panels enlarge a single event of velocity reversal. 

It reveals further that a reversal of is smeared out over 
many (hundred) collisions; in particular, it is not trig- 
gered by hitting a specific obstacle configuration. The 
dynamics may be divided into two states: "running" 
(state A) and "rattling" (state B). A velocity reversal is 
then equivalent to switching from "running" (A) to "rat- 
tling" (B) to "running backwards" (A). We have mea- 
sured the waiting time distribution for the persistence 
of state A and find an exponential distribution with rate 
r~^. Our analysis indicates that the events of velocity 
reversals realize a Poisson process and occur indepen- 
dently and spontaneously. 

Quantitative information on the translational dynam- 
ics is obtained by considering the velocity autocorrela- 
tion function (VACF) Tpcm{t) = (v(t) • v(0)) /v^, see 
Fig. |5] The most prominent feature is a two-step process 
for dense systems: on the microscopic time scale, the 
VACF decays towards a plateau value of 1/2; the termi- 
nal relaxation from the plateau is exponential, providing 
a definition for the zigzag time t^z- Thus, the initial mo- 
mentum has two relaxation channels: rattling in the tube 
and the zigzag motion. 

Let us introduce the correlation functions in the body 
frame := (vait) v„(0)) / (v^) for a =G {||, ±} 

and the orientation correlation Ci{t) = (u(t) • u(0)). 
The decay rates of the three correlation functions are 
well separated at high densities, see inset of Fig. [5] The 
strong suppression of rotational motion in particular al- 
lows to factorize the VACF, 

Am{t) = {[viiit)viiiO)+v±{t)v±{0)] u{t) • uiO))/v^ 

c^l[ij\\{t) + ij±{t)]Ci{t). (3) 

We have checked numerically that this approximation is 
excellent for densities n* > 100. The fastest process is 
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FIG. 5: Two-step relaxation of the CM velocity autocorrelation func- 
tion tpcm{t)- Rescaling of time with the collision rate collapses the 
data for low densities (n* < 2). In the dense regime, a universal 
intermediate plateau develops at 1/2 for Tcoii f <S Tzz; it is pre- 
ceded by a dip reflecting the transversal "rattling in the cage". Inset: 
correlation functions of the different DOFs entering Eq. 

the exchange of transversal momentum with the obsta- 
cles by individual collisions, ip±{t) « exp(— f/rcoii); the 
negative dip after 2 or 3 collisions represents the rattling 
of the needle in its tube. The persistent orientation yields 
the slowest relaxation process, Ci{t) ~ exp(— t/rrot) 
where Trot = l/Dmi. Measuring the rate of longitudi- 
nal relaxation, "011 (*) — exp(— t/rz^), shows an increase 
'Tzz ~ n'^', this provides additional evidence for the en- 
hancement of diffusion upon exploiting a Green-Kubo 
relation, Dcm = (^^^/2) J^Am{t)dt ~ v^Tzz/4:. In 
summary, a ballistic needle between fixed obstacles ex- 
hibits a hierarchy of non-trivial time scales in the dense 
regime, 

Tcoll ~ n"^, Td ~ n°, ~ n'^, r^ot ~ n^, (4) 

where ^ L/v is the disengagement time needed to 
leave the tube [ 1?]. The scaling of all these time scales is 
nicely reproduced by our data; the rapid increase of Trot 
in particular poses the computational challenge to follow 
individual trajectories over many decades in time. 

The interplay of translational and rotational DOFs is 
nicely exhibited in the non-Gaussian parameters (NGP) 
for CM and rotational motion, which are defined as 

aTit) = 1 {(|AR(t)r> / 2 {AR{tf)' - l} , (5a) 

afit) = {A^{t)') I 3 (A^(t)2>' - 1. (5b) 

The persistent momentum renders a^-^ii) negative: after 
equilibration of the transversal momentum, a|"'(t) ex- 
hibits a plateau over several decades for Tcou <C t <C Tzz- 
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FIG. 6: The non-Gaussian parameters for translational and rotational 
DOFs, Eq. (|5}, show a rich structure at high densities including the 
confinement in a tube, the zigzag motion, and translation-rotation 
coupling. Density increases from left to right, and arrows indicate 
the relevant time scales for n* = 800. 



obstacle. This suggests to use ^ t^N above, where 

~ ~ n^^/^. Thus one finds for the exponent 
C = S/A, quite close to the observed ( = 0.8. 

In conclusion, we have observed the conjectured 
power-law divergence of Dcm in a 2D representation of 
needle liquids, thus demonstrating the relevance of our 
results for needle liquids and granular gases. We have 
measured the exponent C, and find its value significantly 
larger than anticipated. The origin of the divergence is 
attributed to the zigzag motion of the needle; further, the 
scaling argument for the exponent C, does not rely on the 
dimension and on the conservation of the single-particle 
energy and thus should hold also for 3D liquids. 
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Noting that the component is uniformly distributed 
for an ergodic sampling of the energy ellipsoid and ap- 
proximating |AR(t)| ~ iiyt, the plateau value is calcu- 
lated to —2/15 in nice agreement with our data. Interest- 
ingly, the plateau is followed by a positive maximum; its 
rising edge scales with Tzz, while the falling edge scales 
with Trot. Such a feature implies that the CM motion is 
non-trivially coupled to the rotational motion. The orien- 
tational confinement of the needle is reflected by a peak 
in a2°'(t) near t ^ TcoU- Then, the signal slowly degrades 
to a shoulder around tj. For longer times, <C t <C Tj-ot, 
the rotational NGP decays slower than l/t, which cor- 
responds to an increasing Burnett coefficient — familiar 
from, e.g., the Lorentz model |22]. 

The value of the exponent C, calls for a theoretical 
explanation. Let the needle initially be oriented along 
the X-axis with most kinetic energy in the longitudinal 
DOF, see Fig. [T] Since the momentum transfer per col- 
lision is perpendicular to the needle, its projection on 
the X-axis is small and fluctuates rapidly around zero, 
Avx ~ sin((/9). Assuming a random walk of Vx, we esti- 
mate the number of collisions N required for a velocity 
reversal to ~ 1/ sm{ip)'^. Accepting that the colli- 
sions occur with rate r~j^j, we have t^z ~ TcohA^. At 
the zigzag scale r^z, the angle of orientation if diffuses, 
and one may approximate ip"^ « 2DmtTzz ^ 1 and thus 
sin(99) « ip. Collecting, we obtain Tzz ~ TcoiiZ-DrotTzz. 
hence ( = 1/2. This reasoning condenses the more so- 
phisticated argument for needle liquids from Ref. 9. Our 
data in Fig. |2l however, show clearly that ( = 1/2 un- 
derestimates the enhancement of diffusion. The above 
argument can be modified by observing that only uncor- 
related collisions contribute to Tzz. The mean free path 
to encounter an obstacle is measured by fTcoii, but the 
needle has to traverse the distance ^ to encounter a new 
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